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1.  Introduction 


1.1  Computation  of  Propagating  Pulses 

The  main  requirement  of  the  Local  Parabolic  Method  (LPM)  development  is  a 
method  to  compute  thin  propagating  wave  equation  pulses  that  do  not  spread  due  to 
numerical  effects. 

(This  part  of  the  research  will,  by  itself  be  important  for  computing  long  distance 
propagation  of  pulses  in  the  time  domain.) 

A  general  method  is  described  to  efficiently  simulate  these  wave  equation  pulses, 
in  Eulerian  computations  on  fixed,  coarse  grids.  The  method,  “Lattice  Confinement”, 
involves  treating  the  features  as  solitary  waves  that  obey  nonlinear,  difference  equations, 
which  are  different  from  Taylor  expansion-  based  discrete  approximations  to  the 
governing,  partial  differential  equations  (pde’s).  These  equations  are  rotationally 
invariant  generalizations  to  multiple  dimensions  of  1-D  discontinuity  confinement 
schemes.  The  method  is  a  generalization  of  an  earlier  method,  “Vorticity  Confinemenf’ 
[1],  which  was  successful  in  efficiently  treating  thin,  vortical  regions. 

For  long  distance  propagation  of  pulses,  direct  discretization  and  solution  of  the 
governing  partial  differential  equations  (pde’s)  using  conventional  Eulerian  Taylor 
expansion  -based  numerical  methods  to  resolve  the  thin  features  can  be  prohibitively 
expensive.  Even  adaptive  unstructured  grid  methods  are  very  expensive  and  complex  for 
general  problems  with  many  small-scale,  time-dependent  features.  Fortunately,  the  details 
of  the  internal  structure  of  these  small  features  are  often  not  as  important  as  integral 
quantities.  The  quantities  of  importance  for  our  purpose  are  the  centroid  motion  and  total 
integrated  amplitude  at  each  point  along  the  pulse  surface.  The  main  issue  in  computing 
these  cases  is  that  conventional  pde-based  methods  require  a  relatively  large  number  of 
grid  cells  (4-8)  across  each  small  dimension  to  treat  a  feature,  such  as  a  wave  equation 
pulse.  Even  then,  the  details  of  the  internal  structure  will  be  mainly  determined  by  the 
discrete  numerics,  and  not  the  physics  of  the  pde.  Also,  numerical  discretization  errors 
will  still  build  up  over  long  distances,  causing,  for  example,  large,  unphysical  spreading. 

This  leads  us  to  the  idea  of  simulating  or  “modeling  “  the  thin  features  directly  on 
the  grid  with  difference  equations,  rather  than  attempting  to  accurately  discretize  the 
pde’s  for  them  using  finite  difference  approximations.  The  idea  of  modeling,  or  solving 
for  small-scale  features  directly  on  the  grid  without  using  smoothness  assumptions  or 
Taylor  expansions,  i.e.,  as  “weak  solutions”,  goes  back  to  work  of  Lax  and  others  [2],  but 
was  applied  mostly  to  shocks.  Shocks,  however,  effectively,  “capture  themselves” 
because  they  have  converging  characteristics.  Harten  [3]  did  treat  contact  discontinuities 
in  this  way,  which  do  not  have  converging  characteristics,  but  for  1-D  compressible  flow. 

First,  the  Lattice  Confinement  method  will  be  reviewed  for  wave  equation  pulses. 
Some  results  will  then  be  presented  for  short  scalar  pulses  obeying  the  wave  equation 
(with  Lattice  Confinement).  Previous  results  [4]  for  a  2-D  pulse  reflecting  from  planar 
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surfaces  will  be  shown,  then  3-D  results  for  a  pulse  reflecting  from  complex  objects  (a 
missile  and  an  aircraft).  This  missile  case  will  involve  multiple,  curvilinear  grids. 

1.2  Main  Features  of  Lattice  Confinement 

(Parts  of  this  section  are  related  to  Ref.  [4]) 

The  main  idea  is  to  treat  thin  features  as  nonlinear  solitaiy  waves  that  “live” 
directly  on  the  grid  lattice,  spread  over  only  a  few  cells.  The  internal  structure  is 
determined  by  the  discrete  lattice  equations.  The  total  amplitude,  centroid  and  (in  a  future 
extension),  a  few  moments,  however,  are  determined  by  the  physics.  These  quantities  are 
transported  across  the  grid  with  essentially  no  numerical  errors  (for  constant  slowness). 
For  smoothly  vaiying  slowness,  small  discretization  effects  will  occur,  based  on  the 
longer  length  scale  of  the  slowness  variation,  but  not  based  on  the  shorter  length  scale  of 
the  pulse.. 

Essentially,  these  discrete  equations  define  a  simple,  implicit  model  which  obeys 
a  “fast”  dynamics,  relaxing  to  an  asymptotic,  propagating  state  in  a  few  time  steps.  In  this 
way  we  can  simulate  the  most  important  physical  effects  of  the  small  scales,  which 
cannot  be  accurately  computed  by  just  discretizing  the  governing  pde’s  on  the  given  grid. 
These  effects  can  be,  for  example,  that  thin  wave  equation  pulses  propagate  over  long 
distances  without  spreading  in  a  smooth,  slowly  vaiying  external  field,  and  that  ftiey  can 
merge  or  reflect,  respectively,  and  thus  change  topology. 

1.2.1  Stationary  Case 

The  formulation  presented  here  is  related  to  that  presented  in  [5]  in  1-D.  First  a 
stationary  pulse  is  discussed,  requiring  only  an  iteration  of  the  Lattice  Confinement 
terms,  so  that  the  simple  asymptotic  form  can  be  seen.  Advection,  in  general,  will  change 
this  form  somewhat.  However,  in  the  limit  of  small  advection  time  step,  or  if  a  number  of 
these  “Confinement”  steps  are  taken  for  each  convection  step,  then  the  following  form 
should  result.  The  same  is  true  for  the  wave  equations.  Results  very  close  to  these  are  also 
found  with  advection  steps  that  are  not  small;  these  are  shown  in  Ref.  [5]. 

We  start  with  an  iteration  sequence  for  a  single-signed  scalar,  0 : 

d,<l)=^-V2\ji(l>-£^\ 

A  t 
or 

<j>n+l  =<pn +h2V\^i<pn -e cD”)  (1.1) 

where  O  is  a  nonlinear  function  of  0  (given  below)  and  n  is  a  diffusion  coefficient  that 
can  include  numerical  effects  in  a  convection  or  wave  equation  solution  (we  assume 
physical  diffusion  is  much  smaller).  The  discretized  grid  cell  size  is  h  and  time  step,  At. 
For  the  last  term,  e  is  a  numerical  coefficient  that,  together  withal,  controls  the  size  and 

time  scales  of  the  confined  features.  For  this  reason,  we  refer  to  the  two  terms  in  the 
brackets  as  “Confinement  terms”. 
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The  two  (positive)  parameters,  e  and  p ,  are  determined  by  the  two  small  scales 
of  the  computation,  h  and  At,  since  we  want  the  small  features  to  relax  to  their  solitary 
wave  shape  in  a  small  number  of  time  steps  and  to  have  a  support  of  a  small  number  of 
grid  cells.  Thus,  even  though  h  may  be  small,  the  Laplacian  will  be  large  and  the  total 
effect  also  large. 

For  the  propagating  pulse  problem,  it  is  assumed  that  the  slowness  field  is 
propagating  is  slowly  varying  compared  to  these  scales  (this  is  required  if  the  grid  cell 
size  and  time  step  are  to  resolve  the  pde’s  governing  this  outer  flow).  Relaxation  of  this 
assumption  will  be  studied  in  the  proposed  effort.  We  then  have  a  two -scale  problem  with 
the  thin  structure  obeying  a  ‘Fast”  dynamics: 

-  eO  »  0 . 


With  propagation  in  a  “slow”,  smooth  external  field,  this  relation  is  still  approximately 
satisfied,  as  verified  by  computations  and  heuristic  arguments  [4]. 

There  are  many  possibilities  for  O .  A  simple  class  is 


O" 


"I  c,»V" 

'  10 

/ 

r  =i<m”  + s . 


(1.2) 


The  above  sum  is  over  a  set  of  grid  nodes  near  and  including  the  node  where  ®  is 
computed.  The  absolute  value  is  taken  and  S  ,  a  small  positive  constant  (~  1CT8 ),  is  added 
to  prevent  problems  due  to  finite  precision.  The  coefficients,  C,,  can  depend  on  /,  but 
good  results  for  many  cases  are  obtained  by  simply  setting  them  as  well  as  pto  1.  Then, 
<l>  is  the  harmonic  mean  of  0  on  the  local  stencil.  Other  forms  could  also  be  used,  with 
p>l.  p  =  oo  corresponds  to  the  minimum  of  the  absolute  value:  for  2-D  and  3-D 

applications  discontinuous  operators  such  as  “min”  will  not  result  in  as  smooth 
distributions  as  continuous  ones,  and  we  use  only  p  =  1  or  p  =  2. 

An  important  feature  of  Lattice  Confinement  is  that  all  terms  are  homogeneous  of 
degree  1  in  Eqn.  1.1  (as  they  are  in  the  wave  equation).  This  is  necessary  because  the 
Lattice  Confinement  terms  should  not  depend  on  the  scale  of  the  quantity  being  confined. 
Another  important  point  is  that  wavelengths  longer  than  the  thin  features  that  are  to  be 
confined  must  have  a  negative  diffusive  behavior,  so  that  the  features  remain  confined, 
that  is  stable  to  perturbations  against  spreading.  This  means  that  d>  must  be  nonlinear:  It 
is  easy  to  show  by  Von  Neumann  analysis  that  a  linear  combination  of  terms,  for  example 
of  second  and  fourth  order,  cannot  lead  to  a  stable  Confinement  for  any  finite  range  of 
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coefficients:  any  wavelength  that  exhibits  negative  diffusion  would  then  eventually 
diverge. 

1.2.2  Wave  Equation  Formulation 

We  start  with  the  1-D  scalar  wave  equation  with  constant  wave  speed,  c,  for 
simplicity.  As  in  scalar  convection,  we  add  an  additional  term  to  control  the  shape  of  a 
short  pulse: 


d2<j)  =  c2320  +  d2xyr 

or,  using  a  simple  time  discretization, 

0"+1  -20"  +0"-'  =c2Atd2x(j)  +  Atd2x\i/  (1.3) 

It  was  seen  in  Ref.  [6]  that  the  addition  of  “Lattice  Confinement”  terms  in  the  form  of 
second  derivatives  of  a  function  that  has  finite  support  do  not  change  the  propagating 
speed  (nor  the  total  amplitude)  of  a  propagating,  confined  pulse.  The  same  is  true  for  the 
wave  equation,  if  an  additional  time  derivative  is  applied.  This  means,  of  course,  that  they 
do  not  change  the  motion  of  the  centroid  of  an  isolated  pulse. 

The  main  constraint  on  the  Confinement  term,  y/ ,  is  that  it  force  an  initial 

isolated,  propagating  compact  pulse  with  a  single  maximum  to  remain  compact  and  not 
develop  any  additional  maxima.  We  use: 

Vn=ltS,r-eS,<y(f)  (1.4) 

In  this  term  O  has  the  form  given  by  Eqn.  (1.2)  in  terms  of  its  argument.  We  have 
defined 


8j"  =r-ri 

Results  will  be  given  in  Section  4  using  this  form. 

An  important  feature  of  the  method  is  that  the  waves  do  not  suffer  a  “phase  shift” 
when  they  pass  through  each  other.  This  is  obvious  for  the  equation  we  want  to  simulate 
-  the  linear  wave  equation.  However,  the  Confinement  term  is  nonlinear.  Such  a  phase 
shift  would  show  up  as  a  kink  in  two  waves  in  2  or  3  dimensions  that  are  passing  through 
each  other,  and  can  be  studied  in  detail  in  1-D.  It  turns  out  that  there  is  no  kink,  to 
plottable  accuracy,  as  can  be  seen  in  the  plotted  scattering  results  in  Sec.  4.1.  Results  for 
the  centroid  trajectories  for  2  pulse  passing  through  each  other  in  1-D  are  presented  in 
Fig.  1.  There,  the  computed  centroids  are  plotted  as  solid  lines  and  the  exact  as  dashed 
(the  periodic  boundary  conditions  can  be  seen  in  the  former).  It  can  be  seen  that  there  is 
no  phase  shift  to  plottable  accuracy.  This  lack  of  nonlinear  interaction  persists,  according 
to  our  study,  in  the  limit  of  small  time  step  (2  orders  of  magnitude  smaller  that  that  of 
Fig.  1),  even  though  0(1 02)  confinement  corrections  were  applied  as  the  pulses  were 
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overlapping.  We  attribute  this  to  the  existence  of  another  conserved  variable,  similar  to 
total  energy.  This  computation  was  done  by  Nick  Lynn  of  University  of  Tennessee  Space 
Institute  (UTSI). 

One  other  important  use  of  Lattice  Confinement  for  the  wave  equation  involves 
cases  with  multiple  grids  with  grid  interfaces.  If  only  the  discretized  wave  equation  is 
used,  with  no  Confinement,  reflections  result  from  small  numerical  errors  at  the  grid 
interfaces,  unless  special  care  is  taken.  Lattice  Confinement  completely  overcomes  this 
problem  and  the  (single)  pulse  propagates  across  the  interface  with  no  reflections. 

Some  of  the  advantages  of  using  direct  difference  equations  instead  of  pde’s  for  a 
related  problem  can  be  found  in  Ref.  [7].  There,  a  very  efficient  formulation  is  derived 
for  the  Helmholtz  equation  by  minimizing  the  L2  norm  of  the  error  of  waves  propagating 
with  fixed  1^1 . 

1.3  Vector  Potential  Formulation 

In  a  recent  promising  development  in  our  wave  equation  research,  we  have 
implemented  a  vector  potential,  A ,  as  the  basic  variable,  instead  of  a  scalar  representing 
the  E  and  B  fields.  If  sufficient  time  is  available,  this  formulation  will  be  investigated  in 
detail. 


Even  though  there  is  no  direct  effect  of  this  alternative  formulation  in  classical 
electromagnetics,  there  is  a  major  advantage  in  the  computation  of  pulse  solutions.  This 
is  because,  for  a  thin  pulse,  E  and  B  only  have  a  small  region  of  support,  but  A  extends 
throughout  the  field.  The  E  and  B  representation  of  the  pulse  is  a  spread  delta  function 
across  the  width,  but  the  A  representation  is  a  step  function.  As  a  result,  it  is  much  easier 
to  capture  the  pulse  by  operating  on  A ,  since  it  is  topologically  “trapped”  by  this  field. 
The  argument  is  exactly  the  same  as  in  Vorticity  Confinement,  where  thin  vortical 
regions  (<a)  are  trapped  in  a  velocity  field,  q,  extending  throughout  space.  (This  is  also 
related  to  the  trapping  of  defects  in  other  field  theories).  The  resulting  confinement  terms 
are  exactly  the  same  in  the  two  cases,  with  A  corresponding  to  q  and  B  corresponding 
to  Co. 


The  formulation  is 

B  =  VxA 

S2A  =  c2V2DA  +  /iStV2DA  +  ed,VDxb 
where  b  is  a  local  harmonic  mean  of  B  at  each  grid  point: 

-i-i 


*-ii 


I— 


N 


where  5,  and  VD  denote  discrete  operators  and  Eqn.  1.2  was  used  with  a  sum  over  N 
neighboring  points,  labeled  l .  In  the  formulation  the  Coulomb  gauge  was  used  so  that  a 
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scalar  potential  does  not  appear  in  the  equation  for  B .  Also,  V  •  A  =  0  is  then  enforced. 
This  is  also  directly  analogous  to  the  incompressibility  condition,  V  •  q  =  0 . 

2.  Objectives 

Develop  “pulse  Capturing”  method  to  generate  approximate  short  wavelength 
solution  which  is  useful  by  itself  for  computing  long  distance  amplitudes  and  which  will 
provide  a  coordinate  system  for  the  locally  parabolic  propagation  method. 

3.  Status  of  effort 

Have  developed  a  pulse  propagation  method  for  tracing  the  path  of  a  short 
wavelength  signal  through  a  medium  with  varying  index  of  refraction  and  over  complex, 
reflecting  terrain.  The  method  is  completely  Eulerian  and  does  not  use  any  markers  that 
can  become  sparse  as  the  signal  spreads  in  the  lateral  direction.  On  the  other  hand,  unlike 
conventional  Eulerian  wave  equation  methods,  the  signal  does  not  suffer  degradation  due 
to  numerical  error  such  as  diffusion,  even  though  the  wavelength  is  of  the  order  of  a  grid 
cell  and  the  signal  can  propagate  over  arbitrarily  long  distances.  During  the  contract 
period,  it  has  been  demonstrated,  both  theoretically  and  numerically,  that  there  are  no 
interaction  effects  when  signals  pass  through  each  other,  even  after  a  large  number  of 
interactions.  This  is  true  of  the  actual  linear  wave  equation  that  is  being  simulated,  but 
had  to  be  shown  for  our  numerical  method,  which  has  a  strong  nonlinear  component. 
Progress  has  also  been  made  in  developing  a  reflecting  formulation  that  can  treat 
complex  terrain.  To  be  feasible  for  realistic  cases,  this  treatment  cannot  use  surface- 
conforming  computational  grids,  but  must  be  able  to  use  a  simple  representation  where 
the  surface  is  "immersed"  in  a  uniform  Cartesian  grid,  with  no  requirements  for  complex 
logic  involving  the  geometry  of  the  "cut"  grid  cells. 

4.  Results 

4.1  Prior  to  4/04 

Lattice  Confinement  was  added  to  the  linear  wave  equation  with  simple,  second 
order  central  discretization.  A  uniform  Cartesian  grid  was  used  with  reflection  on  the 
walls  of  a  square  region,  in  2-D  as  well  as  3-D.  In  certain  cases,  it  is  well  known  that  a 
“tail”  will  develop  behind  a  2-D  wave  front,  while  the  front  remains  sharp  (if  it  is  initially 
sharp).  In  the  2-D  computation  with  Confinement  we  can  keep  a  sharp  pulse  and  suppress 
the  tail.  This  tail  is  smooth  and  could,  if  desired,  be  computed  using  standard  CFD  on  the 
grid.  The  goal,  however,  was  to  show  that  a  pulse  can  propagate  over  long  distances  with 
Lattice  Confinement.  The  same  long  distance  propagation  was  observed  in  3-D  where 
there  were  genuine  pulse  solutions. 

4.1.1  Reflection  from  Planar  Surface 

In  Fig.2  2-D  results  are  presented  for  a  128x128  cell  grid.  It  can  be  seen  that  there 
is  no  perceptible  diffusion  of  the  pulse,  even  after  many  reflections.  These  results  were 
also  presented  in  Ref.  [5]. 
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4.1.2  Scattering  from  Missile 

The  first  computation  involved  a  planar  scalar  wave  impinging  on  the  nose  of  a 
missile.  The  plane  of  the  wave  is  parallel  to  the  longitudinal  axis  of  the  missile.  The 
magnitude  of  the  scalar  is  presented  in  planes  parallel  (i.e.,  the  symmetry  plane)  and 
perpendicular  to  the  longitudinal  axis  of  the  missile,  and  on  the  surface  of  the  missile 
body.  Figures  3  and  4  depict  the  intensity  of  the  planar  wave  before,  during,  and  after 
reflection  from  the  missile  surface  in  the  nose  area  from  the  side  and  from  the  front, 
respectively. 

The  second  computation  involved  the  aft-end  of  the  missile.  Figure  5  depicts  the 
same  sequence,  but  in  the  fin  area,  as  seen  from  the  rear.  The  third  computation  involved 

the  entire  missile.  Figure  6  depicts  a  pulse  reflecting  from  the  nose  area  at45°,  in  the 
symmetry  plane,  as  seen  from  the  side. 

Figure  7  depicts  part  of  the  overlapping  computational  grid.  The  confinement 
method  has  been  generalized  to  curvilinear  grids  for  this  case.  Confinement  also  prevents 
reflections  from  the  grid  interface  regions. 

All  of  the  missile  results  involve  a  pulse  confined  to  about  3  grid  cells  except  in 
the  fine  grid  region  very  near  the  surface.  There,  the  method  reverts  to  standard  CFD. 
Also,  it  can  be  seen  in  Fig.5  that,  unlike  in  geometrical  optics,  there  is  a  diffracted 
component.  This  will  most  likely  represent  a  diffracting  pulse  corresponding  to  the  width 
of  file  computational  one.  We  should  be  able  to  make  corrections  to  the  diffracted 
intensity  to  simulate  pulses  of  other  widths.  This  would  result  in  a  very  general  method 
able  to  treat  diffraction  to  first  order  (in  this  short  pulse  limit). 

4.1.3  Scattering  from  Aircraft 

Next,  results  of  confinement  are  presented  for  scattering  from  an  aircraft  shape. 
Here,  a  uniform  Cartesian  grid  was  used  with  a  special  treatment  of  the  surface  boundary 
conditions.  A  level  set  description  of  the  body  was  used,  which  was  early  derived  from  a 
surface  definition  “STL”  file.  The  contours  of  scalar  intensity  are  presented  in  the  cross 
plane  depicted  in  Fig.  8.  With  these  boundary  conditions,  it  can  be  seen  that  scattering 
from  very  complex  shapes  can  easily  be  computed.  In  Fig.  9,  contours  of  the  scalar 
magnitude  representing  the  electromagnetic  field  are  shown  for  three  different  times. 

The  last  results  involve  the  new  vector  potential  computation.  There,  a  single 
contour  can  be  used  to  represent  each  part  of  the  wave  since  the  potential  extends 
throughout  the  field  and  has  a  definite  range  of  values  within  the  pulse.  Preliminary 
results  are  depicted  in  Fig.  10  for  two  times.  Further  work  is  progressing  on  the  new 
technique. 

4.2  Contract  Period 
4.2.1  Cahn-Hilliard  Equation 
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Based  on  comments  by  Fernando  Reitich  and  Bill  Rider,  a  connection  between 
Cahn- Hilliard  (CH)  -  type  equations  and  our  discrete  convection  equation  (eqn.1.1)  was 
investigated.  This  commonality  includes  the  same  linear  convection  terms  and  the 
appearance  of  a  Laplacian  in  front  of  the  non-linear  term.  This  commonality  is  not 
surprising  since  the  CH  equation  was  originally  proposed  as  a  phenomenological 
description  for  the  dynamics  of  thin  fluid  interfaces  in  otherwise  smooth  flows,  and  our 
equation  describes  thin  pulses  propagating  in  otherwise  smooth  fields.  To  derive  the 
continuum  limit  of  our  equation,  a  parameter  limit  was  derived  such  that  the  pulses  were 
spread  over  an  arbitrarily  large  number  of  grid  cells.  (Though  the  equations  are  similar, 
the  resulting  non-linear  term  is  not  exactly  the  same  as  in  the  commonly  used  form  of  the 
CH  equation).  The  derivation  is  described  in  Fig.  1 1. 

4.2.2  Suppression  of  Nonlinearity  Effects  in  Pulse  Interaction 

As  described  in  Sec.  1.3,  to  represent  intersecting  wave  equation  solutions,  when  pulses 
pass  through  each  other,  there  must  be  no  amplitude  exchange  or  phase  shift.  Otherwise, 
the  actual  wave  equation  being  studied  could  not  be  accurately  simulated,  since  it  is 
linear.  However,  a  nonlinear  term  is  required  in  the  discrete  equation  in  order  to  create  a 
solitaiy  wave  representation  of  the  pulse  which  will  be  non-diflusing.  A  nonlinear  term  in 
the  wave  equation,  of  course,  usually  results  in  interactions  between  intersecting  pulses, 
including  the  above  effects.  As  explained  in  Sec.  4.1,  it  had  been  shown  numerically,  to 
plottable  accuracy,  that  for  our  formulation  these  effects  are  absent.  A  new  analytic  result 
was  obtained  that  shows  that,  for  the  formulation  used,  this  interaction  effect  vanishes  in 
the  Bom  approximation.  This  vanishing  is  due  to  the  fact  that  both  the  Laplacian  and  the 
time  derivative  operator  operate  on  the  nonlinear  term.  This  derivation  is  described  in 
Fig.  12. 

4.2.3  Reflection  from  Complex  Terrain 

Physical  terrain  consists,  of  course,  of  many  irregularities.  For  short  pulses,  or 
high  frequencies,  these  irregularities  can  be  larger  than  the  wavelengths.  This  seems  to 
preclude  attempting  to  achieve  high  accuracy  after  reflection,  except  for  special  cases 
involving  very  smooth  surfaces.  For  this  reason,  one  of  the  main  reasons  for  the  use  of 
higher  order  discretizations  in  conventional  treatments  of  long  distance  propagation  may 
not  involve  increasing  accuracy  after  reflection,  but  to  attempt  to  reduce  numerical 
diffusion  in  the  propagation.  Since  we  do  not  have  this  diffusion  problem,  and  we  are 
treating  complex,  irregular  terrain,  we  feel  that  we  can  use  more  efficient,  lower  order 
discretizations  with  no  loss  of  accuracy.  We  can  then  implement  a  very  effective  method 
for  treating  reflections  (and  absorption).  This  method  does  not  require  complex,  surface- 
fitted  grids,  but  allows  the  terrain  surface  to  be  simply  “immersed”  in  a  uniform  Cartesian 
grid.  This  method  employs  a  “level  set”  representation  of  the  surface  and  can  easily 
accommodate  very  complex  topography  with  little  computational  penalty.  A  picture  of 
the  representation  for  a  simple  object  is  shown  if  Fig.  13.  During  the  computation,  as 
explained  in  Ref.  [4],  each  time  step,  the  wave  amplitude  is  simply  set  to  zero  inside  the 
surface.  The  nonlinear  “confinemenf’  term  keeps  the  surface  definition  sharp.  However, 
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since  the  term  steepens  the  pulse  only  in  the  normal  direction,  the  pulse  surface  is  treated 
as  smooth  in  the  tangential  direction.  As  a  result,  grid  effects  such  as  “staircase”  are 
avoided  and  an  accurate  reflection  is  obtained.  The  scattering  results  presented  in  Sec.  4.1 
used  this  method. 

New  results  have  been  obtained  with  this  method  for  scattering  of  pulses  from  2- 

D  and  3-D  terrain.  The  2-D  results,  including  some  diffraction  effects,  are  presented  in 

Fig.  14.  (Amplitude  contours  are  shown  in  this  and  other  plots  depicting  pulse  scattering). 
Reflection  from  a  3-D  wedge,  depicted  in  Fig.  15,  is  shown  in  the  short  wave  limit  (with 
no  diffraction)  in  Fig.  16. 

4.2.4  Computation  of  Eikonal  Phases  (Frequency  Domain) 

Even  though  the  pulse  representation  in  the  short  pulse  limit  will  be  of  direct  use 

in  computing  amplitudes  and  arrival  times  in  that  limit,  the  computation  of  diffraction 

and  frequency  domain  amplitudes  will  also  be  important. 

The  first  step  towards  this  goal  is  to  use  the  pulse  arrival  time  at  each  grid  node  to 
determine  the  Eikonal  phase,  for  each  pulse  arrival.  Preliminary  results  for  scattering 
from  a  complex  terrain  are  shown  (as  time  contours)  in  Fig.  17  for  first  arrival  times, 
including  a  possible  use  for  computing  diffraction.  In  Fig.  18,  the  direct  diffracted  field  is 
suppressed,  but  both  first  and  second  arrival  times  are  displayed.  Initial  results  for  a 
simpler  case  with  two  sources  and  no  scattering  are  shown  in  Sec.  4.1. 
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11.  Figures 


Figure  1.  Comparison  of  Computed  Centroid  of  Pulses  with  Exact  Solution 
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Figure  5.a~d  Pulse  reflection  from  missile  fins 
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Figure  6.  Pulse  reflection  from  missile 
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Abstract  A  new  method  is  described  that  has  the 
potential  to  greatly  extend  the  range  of  application  of 
current  Eulerian  time  domain  electromagnetic  or  acoustic 
computational  methods  for  certain  problems. 

The  method  involves  adding  a  simple,  nonlinear  term  to 
the  discretized  wave  equation.  As  such,  it  does  not  require 
major  restructuring  of  methods  or  codes  that  have  already 
been  developed.  Researchers  and  engineers  who  are 
solving  problems  for  scattering  or  propagation  of  short 
pulses  should  be  able  to  use  the  new  technique,  in  many 
cases  as  a  simple  “add  on”  or  callable  subroutine,  to  allow 
the  propagation  of  short  pulses  over  long  distances,  even 
if  their  solver  is  low  order  and  the  grid  is  coarse  compared 
to  the  pulse  width  (which  it  must  be  if  the  distances  are 
large).  The  method  has  many  of  the  advantages  of 
Green’s  Function  based  integral  equation  schemes  for 
long  distance  propagation.  However,  unlike  these 
schemes,  since  it  is  an  Eulerian  finite  difference 
technique,  it  allows  short  pulses  to  automatically 
propagate  through  regions  of  varying  index  of  refraction 
and  undergo  multiple  scattering. 

The  new  method,  “  Confinement”  ,  is  based  on  an  earlier, 
very  successful  technique,  “  Vorticity  Confinement”  ,  that 
can  also  be  thought  of  an  “add  on”  ,  which  allows  the 
propagation  of  thin,  concentrated  vortices  over  arbitrarily 
long  distances,  yet  keeps  the  Eulerian  finite  difference 
property  of  the  original  fluid  dynamic  solution  method. 

In  the  paper  the  application  of  Confinement  to  the  scalar 
wave  equation  in  1,  2  &  3  dimensions,  including 
scattering  will  be  described. 

keywords:  Numerical  analysis,  wave  equation, 

computational  acoustics,  computational  electromagnetics. 


1  Introduction 

There  are  many  important  problems  where  thin, 
concentrated  pulses  must  be  numerically  convected  over 
long  distances.  Examples  include  acoustic  and  EM  pulses 
scattered  or  produced  by  aircraft,  rotorcraft  and 
submarines.  Often,  for  these  cases,  the  main  interest  is  in 
the  far  field,  where  the  integrated  amplitude  through  the 
pulse  at  each  point  along  the  pulse  surface  and  the  motion 
of  the  centroid  surface  are  important,  rather  than  the 
details  of  the  internal  structure.  In  general,  these  pulse 
surfaces  can  originate  in  many  places,  multiply  scatter, 
propagate  through  varying  index  of  refraction,  and  have 
complex  topology.  Accordingly,  we  consider  Eulerian 
methods  where  very  general  topologies  can  be  treated. 

Within  this  scope,  there  have  been  many  efforts  over 
decades  to  discretize  and  solve  the  time  dependent  wave 
equations.  Elaborate  codes  have  been  developed  to  treat 
complex  geometries,  such  as  entire  aircraft  (we  have  in 
mind  codes  developed  by  M.  Visbal  of  WPAFB,  V. 
Shankar  of  Hypercomp  Inc.  and  others).  The  application 
of  these  is,  of  course,  limited  by  the  requirement  that 
sufficient  number  of  grid  cells  must  span  the  pulse  to 
accurately  solve  the  equations. 

A  new  method  has  been  developed  that  has  the  potential 
to  greatly  extend  the  range  of  application  of  these 
computational  methods  for  certain  problems.  The  goal  of 
this  effort  is  that  researchers  and  engineers  who  are 
solving  problems  for  scattering  or  propagation  of  pulses 
should  be  able  to  use  the  new  technique,  in  many  cases  as 
a  simple  “add  on”  or  callable  subroutine,  to  allow  the 
propagation  of  short  pulses  over  long  distances,  even  if 
their  solver  is  low  order  and  the  grid  is  coarse  compared 
to  the  pulse  width  (which  it  must  be  if  the  distances  are 
large).  The  new  method  has  many  of  the  advantages  of 
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Green’s  Function  based  integral  equation  schemes  for 
long  distance  propagation.  However,  unlike  these 
schemes,  since  it  is  an  Eulerian  finite  difference 
technique,  it  allows  short  pulses  to  automatically 
propagate  through  regions  of  varying  index  of  refraction 
and  undergo  multiple  scattering. 

The  new  method,  “  Confinement”  ,  is  based  on  an  earlier, 
very  successful  technique,  “  Vorticity  Confinement”  ,  that 
can  also  be  thought  of  as  an  “add  on”  ,  and  which  allows 
the  propagation  of  thin,  concentrated  vortices  over 
arbitrarily  long  distances,  yet  keeps  the  Eulerian  finite 
difference  property  of  the  original  fluid  dynamic  solution 
method. 

Confinement  involves  treating  a  thin  feature,  such  as  a 
pulse,  as  a  type  of  weak  solution  of  the  governing  partial 
differential  equation  (pde).  Within  the  feature,  a  nonlinear 
difference  equation,  as  opposed  to  finite  difference 
equation,  is  solved  that  does  not  necessarily  represent  a 
Taylor  expansion  discretization  of  the  pde.  The  approach 
is  similar  to  shock  capturing  [Lax(1957)],  where 
conservation  laws  are  satisfied,  so  that  integral  quantities 
such  as  total  amplitude  and  centroid  motion  are  accurately 
computed  for  the  feature.  A  more  general  approach  is 
needed,  however,  than  for  shocks,  as  discussed  below. 
Basically,  we  treat  the  features  as  multi-dimensional 
nonlinear  discrete  solitary  waves  that  “  live”  on  the 
computational  lattice.  These  obey  a  “confinement” 
relation  that  is  a  generalization  to  multiple  dimensions  of 
some  earlier  1-D  contact  discontinuity  capturing  schemes. 

Differences  between  Confinement  and  conventional  1-D 
shock  capturing,  are  that: 

First,  unlike  shocks,  characteristics  do  not  point  into  the 
feature,  and  extra  terms  must  be  designed  to  prevent  it 
from  spreading  due  to  numerical  effects  in  the  convection. 
(Harten  [Harten(1978)]  developed  such  a  scheme,  but  for 
contact  discontinuities  in  1-D  compressible  flow.) 

Second,  thin  wave  equation  pulses,  vortex  filaments  or 
thin  streams  of  passive  scalars,  are  intrinsically  multi¬ 
dimensional:  A  concatenation  of  1-D  “  capturing” 
operators  along  separate  axes  will  not,  generally,  give 
smooth  solutions.  Due  to  the  multidimensional  nature,  it 
seems  necessary  to  pay  some  attention  to  the  (modeled) 
structure  within  the  feature,  even  though  it  is  sampled  on 
only  a  few  grid  cells  in  the  cross-section. 

First,  a  short  critique  of  conventional  methods  for  these 
problems  is  given.  The  basic  new  method  is  then 
described.  Initial  results  in  1,  2  &  3D  are  finally 
presented. 

The  method  presented  has  a  similar  goal  to  that  of 
[Bleszynski,  Bleszynski  and  Jaroszewicz  (2004)]  in  that 
they  propagate  a  continuous  wave  surface  in  the  high 
frequency  limit.  The  main  difference  is  that  they  use  a 
system  of  coupled  rays  and  we  use  an  Eulerian  approach. 


Also,  they  are  already  treating  diffractive  effects,  which 
we  are  now  starting  to  do. 

2  Current  Methods 

Conventional  Eulerian  approaches  to  the  wave  equation 
problem  involve,  of  course,  formulating  governing  pde’s, 
discretizing  them  and  solving  them  as  accurately  as 
possible  on  feasible  computational  grids,  assuming 
smooth  enough  solutions.  For  smooth,  non-thin  pulses, 
these  methods  are  well  known  to  converge  to  the  correct 
solution  as  the  number  of  points  across  the  pulse,  N, 
becomes  large:  Error  estimates  are  asymptotic  in  N.  For 
accurate  solutions,  even  higher  order,  complex 
discretization  methods  typically  require  N  to  be  at  least  ~8 
or  1 0  so  that  the  error  obeys  the  large  N  estimate  and  is 
small  [Visbal  and  Gaitonde  (1998)].  Even  then,  solutions 
degrade  over  long  convection  distances  (thousands  of 
pulse  widths).  As  a  result,  conventional  methods  may  be 
inefficient  (or  not  even  feasible)  for  thin  pulses 
convecting  over  long  distances.  Further,  adaptive, 
unstructured  grids  cannot  improve  the  resolution 
significantly  for  realistic  problems  with  many  thin,  time 
dependent  pulse  surfaces. 

3  Confinement  Approach 
3.1  Basic  Features 

For  the  above  reasons,  for  the  problems  considered,  it  is 
important  to  have  only  very  few  (2  or  3)  grid  points  to 
represent  the  cross  section  of  a  pulse  surface  at  each  point 
along  the  surface  and  to  propagate  it  with  no  numerical 
spreading.  This  small  number  of  grid  points  is  consistent 
with  the  desire  to  only  compute  a  few  integral  quantities 
across  the  pulse,  such  as  total  amplitude  and  centroid 
position,  and  perhaps  width  or  a  small  number  of 
moments.  Then,  the  difference  scheme  can,  effectively, 
serve  as  a  simple,  implicit  “solitary  wave”  model  that 
represents  the  wave. 

An  important  point  is  that  both  the  solitary  wave  pulse 
thickness  and  the  physical  pulse  thickness  (they  may  be 
different)  are  assumed  to  be  small  compared  to  the  other 
scales  in  the  region  where  the  method  is  used.  Thus,  the 
pulse  propagates  according  to  geometrical  optics  (high 
frequency  limit)  in  the  region. 

The  basic  idea  is  that  we  want  to  propagate  the  minimum 
amount  of  information  necessary  to  describe  the  pulse. 
When  it  is  thick,  compared  to  other  dimensions,  such  as 
the  nearby  details  of  the  scatterer,  we  may  choose  to  use  a 
fine  grid  and  represent  the  full  physical  pulse  profile.  As 
it  propagates  away,  we  may  just  be  interested,  as 
explained,  in  integral  quantities  at  each  point  along  the 
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pulse  surface,  i.e.,  along  a  ray  normal  to  the  surface.  As 
the  pulse  propagates  away,  we  may  have  to  use  a  coarser 
grid  that  may  even  have  cells  larger  than  the  physical 
pulse  thickness,  while  retaining  this  information  in  our 
“  representative”  solitary  wave. 

An  important  point  is  that,  when  the  pulse  thickness  is 
much  less  than  the  radius  of  curvature  of  the  pulse 
surface,  it  is  more  efficient  to  describe  the  pulse  profile  by 
a  number  of  “  moment  fields”  .  The  resolution  of  the 
thickness  profile  then  depends  linearly  on  the  number  of 
these  moment  fields,  which  only  increases  linearly  with 
the  resolution.  This  is  then  also  true  of  computational 
storage  and  work  requirements.  This  should  be  contrasted 
with  conventional  discrete  Eulerian  schemes,  where  the 
cell  size  is  determined  by  the  required  resolution.  There, 
for  general  configurations  of  surfaces,  the  number  of  grid 
nodes  (and  computer  storage)  in  3-D  increases  like  the 
third  power  of  the  resolution  and,  (including  time  step 
changes),  the  work  increases  like  the  fourth  power. 

As  explained  in  the  next  section,  when  the  grid  is  coarse, 
the  Confinement  method  allows  pulse  surfaces  to 
propagate  over  arbitrarily  long  distances  while  treating 
them  as  nonlinear  solitary  waves,  spread  over  ~  2  grid 
cells,  thus  allowing  information  to  be  accurately 
propagated.  On  the  other  hand,  when  the  grid  is  fine  and 
details  need  to  be  resolved,  the  Confinement  terms 
automatically  become  small  and  the  method  can 
automatically  become  conventional  computational 
acoustics  (or  electromagnetics).  Further,  if  a  pulse 
propagates  through  a  smooth  medium  as  a  solitary  wave 
and  then  encounters  a  scatterer  where  details  must  be 
resolved,  the  pulse  can  be  “  reconstituted”  on  the  (new) 
fine  grid,  if  enough  moments  are  available.  This 
reconstitution  will  require  a  “pulse  shaping”  step.  This 
can  easily  be  effected  since,  in  addition  to  the  common 
positive  numerical  diffusion,  with  confinement,  we  have  a 
stable  total  negative  diffusion,  as  explained  below.  Thus, 
the  fine  grid  pulse  can  be  expanded  or  contracted  until  its 
moments  agree  with  the  correct  values  (this  is  a  subject  of 
current  work).  This  feature  will  be  important  in  many 
cases,  for  example,  when  a  pulse  scatters  from  an  aircraft 
wing,  propagates  many  pulse  widths,  and  scatters  again 
from  a  tail  where  fine  details  must  be  resolved.  Further, 
multiple  scattering  in  inlets  for  short  pulses  should 
provide  an  important  application. 

3.2  Approach 

The  governing  equation  discussed  here  is  the  discretized 
scalar  wave  equation,  with  an  added  Confinement  term: 
(the  approach  als  o  works  for  vectors  or  tensors,  such  as 
Maxwell’s  equations.) 


=  a  2V20  +  — V2  [jii0  -  eO] 

At 

(i) 

where  (j)  is  the  scalar  amplitude,  O  is  the  index  of 
refraction,  jX  is  a  diffusion  coefficient  that  includes 
numerical  effects  (we  assume  physical  diffusion  is  much 
smaller),  and  the  discretized  grid  cell  size  is  h  and  time 
step,  At.  For  the  last  term,  eO ,  £  is  a  numerical 
coefficient  that,  together  with  fl ,  controls  the  thickness 

and  time  scales  of  the  propagating  pulse.  <J>  will  be 
defined  below.  For  this  reason,  we  refer  to  the  two  terms 
in  the  brackets  as  “  Confinement  terms”  .  We  assume 
conventional,  not  necessarily  high  order  discretizations 
are  used  for  the  differential  operators. 

We  have  found  that,  at  least  in  tests  involving  propagation 
through  regions  of  constant  index  refraction,  the  results 
are  similar,  to  plottable  accuracy,  whether  or  not  the  time 
derivative  is  included  on  the  RHS  of  Eq.  (1).  However, 
since  the  time  derivative  enforces  a  relaxation  to  the 
desired  pulse  shape  (as  explained  below),  we  believe  it 
should  be  included  in  general. 

The  basic  idea  is  that  we  want  the  computed  thin  pulses  to 
maintain  their  profile  and  total  amplitude  as  their  centroid 
surfaces  are  propagated  through  the  field.  (We  want  the 
same  for  separate  pulse  fields  representing  moments.)  The 
requirement  that  they  relax  to  their  profile  in  a  small 
number  of  time  steps  and  have  a  support  of  a  small 
number  of  grid  cells  determines  the  two  parameters,  £ 
and  /X  .  Also,  we  assume  that  the  index  of  refraction  field 
in  which  the  pulse  is  propagating  is  slowly  varying  in 
time  and  space  compared  to  these  scales  (this  is  required 
anyway  if  the  grid  cell  size  and  time  step  are  to  resolve 
this  field).  We  then  have  a  two-scale  problem  with  the 
thin  pulse  obeying  a  “  fasf ’  dynamics. 

V2(/f0-£< I>)  =0, 


Thin  pulses  are  then  propagated  through  the  field  by  the 
“slow”  variable,  (7  .  Exactly  the  same  type  of  discussion 
applies  to  the  convection  of  passive  scalars,  as  described 
in  Ref.  [Steinhoff,  Fan,  Wang  and  Dietz  (2003)]. 

In  general,  the  integrals  that  we  are  interested  in  are  not 
sensitive  to  the  parameters  £  and  [l  over  a  wide  range 

of  values,  as  long  as  the  computed  pulses  are  thin. 

An  important  feature  of  the  Confinement  method  is  that, 
since  it  is  a  second  derivative  in  space  and  first  in  time, 
the  total  amplitude  and  centroid  of  the  surface  are  not 
changed  by  the  added  confinement  terms,  even  under 
discretization. 
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3.3  Formulation 

The  formulation  for  Confinement  will  first  be  described 
for  a  stationary  pulse,  for  clarity.  The  scalar  formulation 
presented  here  is  related  to  that  presented  in  [Steinhoff, 
Wenren,  Underhill  and  Puskas  (1995),  Steinhoff,  Puskas, 
Babu,  Wenren  and  Underhill  (1997)]  in  1-D.  This  “  fast’ 
dynamics  will  be  realized  in  a  wave  equation  computation 
in  the  limit  of  small  time  step,  or  if  a  separate 
“  Confinement”  iteration  is  done  each  time  step.  Excellent 
results  are  found  with  convection  and  are  shown  in 
[Steinhoff,  Fan,  Wang  and  Dietz  (2003)]  for  vorticity  as 
well  as  convecting  passive  scalars. 

For  this  case,  we  have  an  iteration  for  a  non-negative 
scalar,  (j) : 

<pn+x  =<j)n  +iih2V2(pn -eh2V2<i>n 
(3) 


where  the  number  of  terms  h  the  sum  is  N=9.  Here,  we 
assume  0”  >  0 .  Negative  values  can  also  be 
accommodated  with  a  small  extension.  Both  /I  and  £ 
are  positive. 

An  important  feature  is  that  all  terms  are  homogeneous  of 
degree  1  in  Eq.  (4).  This  is  important  because  the 
confinement  should  not  depend  on  the  scale  of  the 
quantity  being  confined.  Another  important  feature  is  the 
nonlinearity.  It  is  easy  to  show  that  a  linear  combination 
of  terms,  for  example  of  second  and  fourth  order,  cannot 
lead  to  a  stable  confinement  for  any  finite  range  of 
coefficients. 

For  smooth  (f)  fields  (long  wavelengths),  the  last  term  in 
Eq.  (1)  represents  a  diffusion.  If  jj.  <  £ ,  the  total 

diffusion  (in  the  long  wavelength  limit)  is  negative. 
However,  the  iteration  of  Eq.  (3)  is  still  stable  and  has 
been  observed  to  converge  for  values  of  £  several  times 
that  of  jl  (depending  on  value  of  fl ). 


where 


(4) 

r=\<t>\n  +8 

(5) 


where  the  sum  is  over  a  set  of  grid  nodes  near  and 
including  the  node  where  is  computed,  the  absolute 
value  is  taken  and  S  ,  a  small  positive  constant  (~  10  8 ) 
is  added  to  prevent  problems  due  to  finite  precision.  The 
coefficients,  Ct ,  can  depend  on  /,  but  good  results  are 

obtained  by  simply  setting  them  all  to  1  for  the  wave 
equation  (different  values  are  used  for  passive  scalar 
convection  to  avoid  using  downwind  values  [Steinhoff, 
Fan,  Wang  and  Dietz  (2003)].  Eq.  (4)  is  related  to  the 
harmonic  mean. 

For  example,  in  2D,  except  for  convecting  scalars,  the 
form  used  in  this  study  is 


3.4  A  nalysis  of  Small  Time  Step  Form 

(Sections  3.4  and  3.5  are  close  to  part  of  Ref.  [Steinhoff, 
Fan,  Wang  and  Dietz  (2003)]. 

Stability  of  the  iteration  as  n  — »  °°  can  easily  be  shown 
[Steinhoff  and  Lynn  (2002)]  for  a  range  of  values  of  /I 
and  £ ,  including  jl<£ .  We  only  have  to  start  with  a 
non-negative  initial  (0°)  field  and  show  that,  for  the  fl 
and  £  values,  (f>n  remains  non-negative.  Since  the  sum 
of  <j)  values  is  conserved,  there  is  thus  an  upper  bound. 
Assuming  convergence  as  «—><»,  we  have 

V2(H<j)-& X>)  =  0 
(7) 

If  <j)  (and  hence  <J> )  vanishes  in  the  far  field,  away  from 
the  pulse,  we  have  [A(f)  =  £& , 

If  the  point  (i,  j)  is  given  the  label  /  =  0 ,  we  then  have 


~2<t>'A=o 

£N  i 


There  are  many  solutions  of  this  equation.  The  ones  of 
importance  to  us  are  of  the  form 


(pij  =  A  sec  h[a(z  -  z0)] 

(9) 

z  =  xt  cos  6  +y  j  sin  6 

(10) 


(6) 
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A,  CC ,  Z0,  6  constant,  and  where  Xt  =  ih  ,  y y.  =  jh ,  h 
is  the  grid  cell  size,  and  we  use  the  form  corresponding  to 
C,=  1  in  Eq.  (4).  This  converges  to  a  straight  pulse  (in 

2-D)  concentrated  about  a  line  at  angle  6 .  It  is  easy  to 
see  that  CC  satisfies 

—  =  [1  +  2  ch(ah  cos  0)  +  2ch(ah  sin  0)]  /  5 
P 

(H) 

for  N  =  5  . 

An  important  point  is  that  we  obtain  close  to  the  same 
invariance  properties  as  the  original  pde:  The  solution  is 
translationally  invariant  (z0  is  arbitrary)  and  close  to 
rotationally  invariant  (0  is  arbitrary  with  a  width,  given  by 
a  in  Eq.(l  1),  having  some  dependence  on  0). 

3.5  Convection  of  Passive  Scalar 

Since  the  wave  equation  is,  of  course,  closely  related  to 
the  convection  equation,  we  present  some  analysis  for  the 
latter,  since  it  is  simpler.  This  analysis  shows  that  the 
pulse  convects  with  the  weighted  mean  velocity,  where 
the  pulse  amplitude  is  the  weight.  This  “  Ehrenfest”  type 
of  relation  should  extend  to  the  full  wave  equation. 

The  following  argument  assumes,  for  each  convection 
step  (n),  there  is  at  least  one  confinement  step  so  that  the 
feature  remains  compact.  If  (J>  represents  a  confined 
passive  scalar,  then,  using  a  conservative  convection 
routine,  we  have  the  following  relationships  for  the 
dynamics  of  the  convecting  solitary  wave  (we  describe 
the  2-D  case  for  simplicity): 

We  have  a  discretization  of 

3,0  =  -V-(§0)  +  /t2V2(/x0-eO)/Af 

(12) 

assuming  V  •  q  =  0  .  Then, 

0M+1  =0"  -A tVdisc  •(#)+^2V^c(^0-eO) 

(13) 

where  discrete  operators  are  labeled. 

For  conservative  discretization,  the  total  amplitude 

n = IX 

'j 

(14) 

is  independent  of  n.  If  we  define  the  centroid 


V 

(15) 

and  the  weighted  mean  velocity 

<?>”-! 

ij 

(16) 

where  is  the  (fixed)  position  vector  of  node  (i,  j),  and 
0^  and  qtj  are  the  scalar  value  and  the  velocity  at  that 
node,  then  the  centroid  evolves  according  to: 

<X>n+1=<X>"  +At<q>" 

(17) 

Since  we  are,  at  this  point,  only  interested  in  the 
“expectation  values”  for  thin  pulses  and  that  the  pulses 
remain  compact,  spread  over  only  a  few  cells,  this 
Ehrenfest-type  relation  is  exactly  what  we  need.  Only  the 
variables  of  importance  are,  effectively,  solved  for.  This 
shows  that  the  pulses,  when  isolated,  evolve  as  surfaces 
with  essentially  no  internal  dynamics  (assuming  they 
remain  confined  as  thin  surfaces).  However,  we  keep  the 
very  important  Eulerian  feature  that  the  number  of  pulses 
is  not  fixed.  We  could,  for  example,  create  additional 
solitary  waves  by  inserting  a  source:  No  additional 
computational  markers  need  be  created,  as  in  Lagrangian 
schemes.  For  this  study,  we  show  that  pulses,  for 
example,  reflect  and  thereby  increase  automatically  in 
number.  This  will  be  seen  in  the  results  of  Sec.  4. 

4  Results 

For  the  scalar  wave  equation,  a  simple  second-order 
centered  difference  method  was  used  for  the  discretization 
of  Eq.  (1).  We  solve  it  through  two  steps:  the  first  is  a 
conventional  wave  equation  solver  step,  and  the  second  is 
the  confinement  step 

0*  =  20”  -0”-'+ctWvL.0" 

(18) 

0”+1  =0*  +/t2V^c5“(/x0*-eO*)  (19) 

where  V^iJC  and  8~  are  the  second-order  centered 

difference  approximation  for  the  Laplace  V2  operator 
and  the  backward  difference  operator  in  n. 
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For  the  1-D  and  2-D  result  plots,  axes  are  labeled  with 
grid  node  location.  For  the  2-D  and  3-D  results,  plots  of 
amplitude  are  made  using  dense  contours. 

4.1 1-D  Wave  Equation 

We  start  a  single  pulse  in  the  center  of  a  256  cell  grid  with 
periodic  boundary  conditions.  This  pulse  has  an  initial 
amplitude  of  2  at  the  single,  central  grid  point  and  is  zero 
at  all  other  points.  The  resulting  two  opposite  moving 
waves  are  shown  in  Fig.  1  for  no  Confinement  (la)  and 
Confinement  (lb).  A  minimal  diffusion  necessary  for 
stability,  jl  =  0.025 ,  with  £  =  0  was  used  in  la  while 

in  lb  H  =  0.1 ,  £  =0.8.  The  rapid  diffusion  can  be 
seen  in  la,  while  in  lb  the  bulk  of  the  pulses  remain 
confined  to  ~3  grid  cells  and  have  no  diffusion.  Other 
tests  show  that  they  continue  unaltered  for  up  to  about 
106time  steps,  and  beyond  if  double  precision  is  used. 
This  was  also  shown  with  an  earlier  Confinement  form  in 
Ref.  [Steinhoff,  Wenren,  Underhill  and  Puskas  (1995)]. 

It  should  be  noted  that  even  though  the  Confinement  is 
nonlinear,  there  is  virtually  no  interaction  when  the  waves 
pass  through  each  other.  This  was  shown  in  detail  in  Ref. 
[Steinhoff,  Wenren,  Underhill  and  Puskas  (1995)]. 


Figure  1:  ID  pulse  propagation 


4.2  2-D  Wave  Propagation 

Waves  were  propagated  on  a  (128) 2  cell  grid  with 
reflecting  boundary  conditions.  Confinement  values  used 
were  fl  =  0.08,  £  =0.6.  Of  course,  the  actual  wave 
equation  exhibits  a  “tail”  behind  a  pulse  in  2-D.  This  can 
be  seen  to  be  suppressed  by  the  Confinement,  and, 
effectively,  only  the  steep  pulse  front  is  accurately 
computed.  The  tail,  since  it  is  smooth,  could  be  computed 
with  no  Confinement.  The  main  interest,  however,  is  in  3- 
D  and  this  was  not  done. 

4.2.1  Convex  Wave 

An  outward  propagating,  initially  circular  pulse  surface 
(diameter  64  cells)  was  computed.  It  can  be  seen  in  Fig. 2 
that  it  remains  sharply  confined,  even  after  many 
reflections.  Again,  as  in  1-D,  there  is  no  discernable 
interaction  between  intersecting  waves. 

4.2.1  Concave  and  Convex  Waves 

The  same  computation  was  done  as  in  4.2.1,  but  with  an 
initially  2:1  elliptical  surface,  with  64  cell  major  axis. 
Both  inward  and  outward  moving  waves  were  formed. 
The  inward  moving  wave  can  be  seen  to  form  cusps  and 
“swallowtails”  .  These  are  only  initially  resolved,  because 
of  the  coarseness  of  the  grid.  The  basic  discrete  wave 
equation  method,  without  Confinement,  would,  of  course, 
not  do  better.  It  is  well  known  that  refinement  is  needed  in 
such  regions  for  direct  application  of  finite  difference 
schemes  [Benamou  and  Solliec  (2000)]. 


4.3  3-D  Convex  Wave  Propagation 


“k" 


Figure  3:  2D  convex  and  concave  wave  propagation: 
Cusp  Formulation 
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Conclusion 
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A  new  Eulerian  technique  is  introduced  for  solving  the 
wave  equation  for  short  pulses.  The  method, 
“  Confinement”  ,  reduces  to  standard  Eulerian  ones  for 
smooth,  long  wavelength  pulses.  However,  unlike 
conventional  schemes,  it  does  not  diffuse  short  pulses. 
Instead,  they  are  “  Confined”  and  propagate  as  nonlinear 
solitary  waves  that  “  live”  on  the  computational  lattice.  As 
such,  they  can  be  propagated  over  indefinitely  long 
distances,  while  remaining  only  2~3  cells  thick.  These 
pulses  represent  the  short  physical  pulses  and  accurately 
propagate  integral  quantities  at  each  point  along  the  pulse 
surface,  such  as  total  amplitude,  centroid  position,  pulse 
width  and  other  desired  moments.  It  is  argued  that,  for 
thin  pulses,  the  method  can  easily  be  implemented  in 
existing  codes,  allowing  them  to  be  extended  to  treat 
much  higher  wavelengths/shorter  pulses,  without 
extensive  reformulation.  Examples  are  shown  in  1-D,  2-D 
and  3-D. 
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Figure  4:  3D  convex  wave  propagation 
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